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Abstract.  Time  accurate  numerical  simulations  of  a  high  temperature  source  cloud  of  gas  expanding  into  an 
ambient  atmosphere  are  performed  using  a  multiple  temperature  gas  model  and  the  direct  simulation  Monte  Carlo 
(DSMC)  method.  The  multi-temperature  approach  uses  continuum  conservation  equations  derived  from  the 
Boltzmann  equation  via  first-order  Chapman-Enskog  expansion  and  zero-order  non-isotropic  velocity  distribution 
function.  These  equations  are  solved  numerically  using  a  kinetic  flux  splitting  method  for  inviscid  fluxes  and  a 
central  difference  scheme  for  the  viscous  fluxes  in  a  time  accurate  manner.  The  DSMC  technique  is  a  well-established 
particle  approach  for  rarefied  flows.  The  source  and  ambient  densities  and  temperatures  are  varied  to  establish  the 
range  of  applicability  and  relative  accuracy  of  the  two  methods.  For  most  of  the  source  and  the  ambient  parameter 
considered,  the  two  methods  give  similar  predictions  of  the  basic  flowfield  evolution,  and  of  the  degree  of 
translational  non-equilibrium  that  arises. 


INTRODUCTION 

The  expansion  of  a  highly  energetic  gas  cloud  into  a  rarefied  environment  is  a  phenomenon  that 
arises  in  many  practical  problems,  including  supernova  explosions,  high  altitude  rocket  plumes,  and 
active  space  experiments  where  interaction  of  the  artificial  plasma  and  neutral  gas  cloud  with  the 
ambient  atmosphere  is  studied.  During  the  rapid  expansion  of  the  gas  cloud,  the  spherical  gas  front  or 
contact  surface  propagates  outward  through  the  ambient  atmosphere,  while  a  rarefaction  wave 
propagates  inward  into  the  source  cloud.  Recent  time-accurate  continuum  simulations1  of  a  high- 
temperature  metallic  gas  cloud  at  high  altitude  with  radiation  transport  show  the  flow  field  generally 
includes  continuum,  transitional,  and  rarefied  regions  depending  on  the  absolute  and  relative  densities 
and  thermal  energies  of  the  source  cloud  and  ambient  atmosphere.  Thus  it  seems  that  particle 
simulation  methods  such  as  the  direct  simulation  Monte  Carlo  (DSMC)  method2  would  be  more 
suitable  for  simulating  these  flows.  However,  the  unsteady  nature  of  these  flows  makes  the  particle 
approach  a  formidable  computational  task.  Other  features  that  add  further  complexity  to  the  particle 
treatment  are  the  presence  of  very  dense  gas  in  the  source  region  and  the  potential  for  phase  change 
as  well  as  coupled  radiative  transport.  Moreover,  the  dilute  gas  assumptions  made  in  these  particle 
simulation  methods  make  them  unsuitable  for  dense  regions  of  the  cloud.  This  suggests  that  at  present 
a  continuum  method,  which  to  some  extent  accounts  for  the  non-equilibrium  nature  of  the  flow,  may 
provide  efficient  and  accurate  simulations  of  the  gas  cloud  expansion. 

The  purpose  of  this  study  is  to  improve  the  fundamental  understanding  of  the  initial  evolution 
of  the  gas  cloud  expansion,  including  the  formation  of  the  shock  wave  in  the  ambient  atmosphere  and 
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the  degree  of  translational  non-equilibrium  that  arises.  Time  accurate  numerical  simulations  of  a 
rapidly  expanding  high  temperature  gas  cloud  are  performed  using  both  a  multi-temperature  gas 
dynamics  model3'4  and  the  DSMC  method2.  The  influence  of  the  source  radius  and  source  and  ambient 
density  are  evaluated.  Comparison  of  results  for  the  two  methods  establishes  the  range  of 
applicability  and  the  relative  accuracy  of  each. 


COMPUTATIONAL  APPROACHES 

Conservation  equations  of  the  multi-temperature  gas  dynamic  model  in  spherical  symmetric 
coordinates  (Ref.l)  are  solved  numerically  in  a  time  accurate  manner.  The  inviscid  part  of  the  flux  is 
calculated  at  each  cell  interface  using  a  kinetic  flux  splitting  method. 1  A  local  ellipsoidal  distribution 
function  of  the  form 
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is  then  defined  for  each  cell  and  used  to  formulate  analytical  expressions  for  the  inviscid  fluxes  in  the 
radial  (r)  and  perpendicular  (p)  directions.  The  c’„  c’p,  and  c’t  are  the  components  of  the  random 
molecular  velocity  in  the  above  expressions.  As  with  other  flux-splitting  techniques,  this  inviscid  flux 
is  formally  first-order  accurate  in  space.  The  second-order  accuracy  in  space  is  achieved  by  using  a 
MUSCL  variable-extrapolation  scheme5  in  conjunction  with  a  minmod  slope  limiter  to  ensure 
monotonicity.  The  viscous  fluxes  are  calculated  with  a  second-order  accurate  central  difference 
scheme.  Fourth-order  accuracy  in  time  is  obtained  using  a  Runge-Kutta  time  integration  method6.  The 
time  step  chosen  is  the  minimum  of  a  Courant-based  time  step  and  a  step  based  on  the  magnitude  of 
energy  exchange  between  the  radial  and  perpendicular  directions. 

The  DSMC  calculations  are  performed  using  a  wholly  conventional  one-dimensional  method  well 
described  in  reference  2.  Some  features  of  the  current  implementation  are  given  in,  e.g.,  reference  7. 
Profiles  of  macroscopic  properties  at  a  given  time  are  generated  as  instantaneous  snapshots  of 
sampled  molecular  quantities. 


RESULTS  AND  DISCUSSION 

The  source  cloud  is  assumed  to  be  spherical  with  uniform  initial  density  and  temperature.  The 
source  is  pure  argon  with  an  assumed  Maxwellian  (inverse  fourth  power  law)  inter-atomic  potential 
(  Tref  =213Kand\iref  =  2.517.tl0~5  kg /ms).  Two  values  of  source  gas  density  p,  at  a  temperature  of  Ts 
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=  3000  K  are  considered  in  the  present  calculations.  The  outer  boundary  condition  is  assumed  to  be 
the  standard  atmosphere  at  the  given  altitude  (70  or  85  km).  To  focus  on  cloud  expansion  and  not  on 
molecular  diffusion  and  other  multi-species  issues,  the  ambient  atmosphere  is  simulated  as  pure  argon 
at  the  ambient  atmospheric  density  and  temperature. 


In  the  multi-temperature  continuum  calculations  the  initial  spherical  source  volume  and  flow 
domain  are  divided  into  uniformly  distributed  spherical  cells  of  constant  width  0.0002m.  One  cell  is 
added  to  the  outer  boundary  after  each  time  step,  allowing  the  grid  to  grow  in  time  as  the  gas  cloud 
expands.  This  procedure  reduces  the  computational  task  during  early  phases  of  the  solution  because 
computation  may  need  very  large  number  of  cells  at  later  times. 

The  DSMC  calculations  use  a  uniform  radial  grid  with  cell  dimension  comparable  to  the  initial 
source  mean  free  path.  Collision  sub  cells  were  utilized  to  further  improve  the  spatial  fidelity. 
Maxwell  molecule  collision  interaction  was  approximated  using  the  Maxwell  limit  of  the  VHS  model. 
Sensitivity  studies  (not  shown)  confirm  the  results  to  be  grid  and  timestep  independent.  Since  spatial 
weight  factors  are  not  used,  initial  molecular  population  of  a  cell  is  controlled  only  by  the  local  gas 
number  density  and  the  cell  volume.  Control  of  statistical  noise  present  in  the  sampled  snapshots  is 
possible  in  this  case  only  by  recourse  to  large  numbers  of  simulated  molecules.  The  typical  simulation 
required  40  million  molecules,  1000  cells,  10  subcells  per  cell,  and  2500  timesteps. 


Figures  1-4  show  the  evolution  of  the  density,  temperature,  and  velocity  profiles  for  a  source 
cloud  expanding  into  a  70  km  altitude  ambient.  The  initial  source  density  and  radius  for  this  case  are 
p,  =  0.005  kg/m3  and  rs  =0.1  m,  corresponding  to  a  source  Knudsen  number  of  Kn,  ~  5.6xl04.  The 
source  to  ambient  density  ratio  is  approximately  60.  The  initial  source  energy  is  19.6  J  with  energy- 
density  of  4.68xl03  J/m3.  These  figures  show  that  a  well-defined  shock  wave  propagates  outward  in 
the  ambient  atmosphere.  The  shock  is  followed  by  a  contact  surface,  which  separates  the  expanding 
gas  from  the  ambient  environment.  A  strong  rarefaction  wave  travels  behind  the  contact  surface 
towards  the  center  of  the  source.  The  final  time  shown  (t=lxl04  s)  corresponds  to  approximately 
one  acoustic  time  (based  on  the  initial  source  radius  and  speed  of  sound). 

Translational  non-equilibrium  arises  in  the  shock  wave  as  the  gas  cloud  expands  (Fig.  3).  The 
ratio  of  the  radial  temperature  to  the  parallel  temperature  (Tr/Ts)  as  a  function  of  normalized  gas 
cloud  radius  is  shown  in  Fig.  3  for  three  different  times.  The  maximum  value  of  this  ratio  is  more  than 
two  which  occurs  in  the  shock  wave  region  (Fig.  3).  This  suggests  that  the  shock  wave  region  is  in  a 
highly  non-equilibrium  thermal  state.  The  majority  of  the  expanding  gas  cloud  remains  in 
translational  equilibrium,  due  to  the  relatively  high  source  density.  A  slight  overshoot  is  observed  in 
the  radial  temperature  just  behind  the  shock  wave  (Fig.  2),  but  no  overshoot  is  present  in  the 
perpendicular  temperature  (not  shown). 

The  radial  velocity  increases  linearly  with  radius  through  the  cloud  for  most  of  the  flow  field 
(Fig.  4).  Also,  the  radial  velocity  gradient  decreases  with  time  as  the  gas  cloud  expands  into  the 
ambient  atmosphere. 
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Altitude  70  km 


Altitude  70  km 


Fig.  1  Density  Profiles  Fig.  2  Radial  Temperature  Profiles 


Altitude  70  km 


Fig.  3  Parallel  Temperature  Profiles 


Fig.  4  Radial  Velocity  Profiles 
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DSMC  calculations  predict  a  slightly  thicker  shock  than  predicted  by  the  multi-temperature 
continuum  model  (Figs.  1-4).  The  DSMC  results  are  expected  to  be  more  accurate  than  the 
continuum  method  in  the  highly  non-equilibrium  leading  edge  of  the  shock.  Also,  the  difference  in 
shock  wave  thickness  is  slightly  more  pronounced  in  the  temperature  profiles  because  temperature  is 
calculated  based  on  higher  moments  of  the  distribution  function.  The  DSMC  results  show  relatively 
large  scatter  near  the  origin  due  to  the  small  cell  volumes  present,  but  this  does  not  obscure 
interpretation  of  the  remainder  of  the  flowfield.  The  overall  qualitative  and  quantitative  comparison 
of  the  flow  field  predictions  by  the  two  methods  is  very  good  for  this  case. 

A  second  set  of  calculations  (Figs.  5-8)  is  also  performed  at  70  km  altitude  for  expansion  of  a 
gas  cloud  with  five  times  smaller  source  density  and  half  the  radius  of  the  previous  case.  For  this  case. 
rs  =  0.05  m  and  ps  =  0.001  kg/m3,  resulting  in  Kns  ~  5.6xl0"3  and  p,  /p„  ~  12.  The  source  energy  is 
0.489  J.  and  the  energy  density  is  934  J/m3.  These  calculations  illustrate  the  effect  of  initial  source 
energy  on  the  gas  cloud  expansion  for  given  ambient  conditions.  Figures  5-8  show  the  same  general 
features  of  the  gas  cloud  expansion  as  are  observed  in  the  first  set  of  calculations  (Fig.  1-4)  except 
that  a  shock  wave  does  not  form  in  the  ambient  because  there  is  not  enough  source  energy  to  raise 
the  expansion  front  velocity  to  high  supersonic  value  (Fig.  8).  Comparison  of  Fig.  6  and  Fig.  7  of  the 
temperature  profiles  shows  considerable  thermodynamic  non-equilibrium  in  the  expansion  front 
region. 


Altitude  70  km 


Altitude  70  km 


Fig.  6  Radial  Temperature  Profiles 
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Fig.  7  Parallel  Temperature  Profiles 


Fig.  8  Radial  Velocity  Profiles 


In  order  to  see  the  effect  of  ambient  rarefaction  on  the  gas  cloud  expansion,  a  third  set  of 
calculations  is  performed  for  the  source  conditions  used  in  the  first  set  (Fig.  1-4),  but  with  an  ambient 
corresponding  to  an  altitude  of  approximately  85  km.  The  ambient  values  of  density  and  temperature 
are  p„  =  8.2196  x  lCT6  kg/m3  and  Ta  =188.893  K,  giving  a  density  ratio  p,  /p„  -  600.  Figures  9  and  10 
show  normalized  plots  of  density  and  radial  temperature  as  a  function  of  normalized  radial  distance. 
Fig.  9  shows  that  the  multi-temperature  model  predicts  a  well-defined  thin  shock  in  the  ambient, 
whereas  the  DSMC  calculations  predict  a  thick  and  diffused  shock  wvae.  The  multi-temperature  model 
cannot  resolve  the  large  non-equilibrium  effects  present  in  the  leading  edge  of  the  disturbance,  and  for 
this  reason  the  radial  temperature  in  the  shock  region  predicted  by  the  multi-temp  model  is  lower 
than  predicted  by  the  DSMC  method  (Fig.  10).  The  difference  in  parallel  temperatures  in  the  shock 
wave  region  (not  shown)  is  much  smaller  than  for  radial  temperatures.  Also,  the  difference  in  radial 
velocity  predicted  by  the  two  methods  in  the  shock  wave  region  (not  shown)  is  small.  The  overall 
comparison  of  flow  field  predictions  by  two  methods  is  reasonably  good. 


Altitude  85  km 


Altitude  85  km 


Fig.  10  Radial  Temperature  Profiles 


CONCLUSIONS 

Computations  carried  out  for  a  variety  of  source  and  ambient  conditions  demonstrate  the  effects 
of  source  energy  and  ambient  rarefaction  on  the  evolution  of  the  gas  cloud.  For  the  low  source  energy 
conditions  at  70  km  altitude,  no  shock  wave  propagates  in  the  ambient  because  the  gas  front  velocity 
is  barely  supersonic.  Higher  source  energy  conditions  at  the  same  altitude  result  in  a  well  defined  shock 
wave,  with  considerable  translational  non-equilibrium  present  in  the  shock  wave  region.  The  multi- 
temperature  continuum  model  and  the  DSMC  method  predict  a  radial  temperature  overshoot  in  the 
gas  front  region  for  the  dense  and  high  energy  source  conditions  at  70  km  altitude.  No  radial 
temperature  overshoot  is  predicted  for  the  same  source  expanding  at  85km  altitude  where  ambient 
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atmosphere  is  more  rarefied  by  an  order  of  magnitude.  The  comparison  of  flow  structure  of  the  gas 
cloud  expansion  predicted  by  the  multi-temperature  model  and  the  DSMC  method  shows  the  multi- 
temperature  model  is  a  reasonably  good  computational  tool  for  low  to  moderate  rarefied  and  non¬ 
equilibrium  conditions. 
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